SERIES REPRESENTATIONS AND SIMULATION OF ISOTROPIC RANDOM FIELDS IN THE EUCLIDEAN SPACE

This document represents the Pytho codes associated with Series Representations and Simulations of Isotropic Random Fields in the Euclidean Space by Zhengwei Ma and Chunsheng Ma

# python imports

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import itertools

from numpy.random import normal, uniform, gamma
from scipy.special import jv
import plotly.graph_objects as go
import multiprocessing as mp

Setup

Here we set some initial values for simulation, as well as create some helper functions

r = 10  
theta = 0.2

simulation_number = 1000 # number of times to simulate

m = 100 # number of points to simulate

r_list = np.linspace(0, 10, m)
theta_list = np.linspace(0, np.pi * 2, m)
R, T = np.meshgrid(r_list, theta_list)  # R, T are matrices of size 100 by 100
X, Y = R*np.cos(T), R*np.sin(T) # express in polar coordinates

def generate_z(function):
    z_list = []
    with mp.Pool(mp.cpu_count()) as p:
        z_list = p.starmap(function, itertools.product(range(m), range(m)))
    
    Z = np.reshape(z_list, (m, m))
    return Z

(i) Example 3.3 (continued), Case I: \nu=0

We want \begin{equation} y_{1} = -log(U_{0, 1}) \end{equation}

\begin{equation} y_{2} = U_{log(2), log(8)} \end{equation}\begin{equation} v_{1} = y_{1}^{0.5} * exp(\frac{-y_{2}}{2}) \end{equation}\begin{equation} w_{1} = U_{0, 1} \end{equation}\begin{equation} u_{1} = U_{0, 1} \end{equation}

def example_1(i, j):
    """
    This function takes in two positional arguments, i and j, which represents the corresponding 
    values in a m by m matrix (based on the value set in cell 4)
    """
    r = R[i, j]
    theta = T[i, j]
    
    z_output = []
    for i in range(simulation_number):
        y1 = -np.log(uniform(low=0, high=1, size=1))
        y2 = uniform(low=np.log(2), high=np.log(8), size=1)
        v1 = (y1 ** 0.5) * np.exp(-y2/2)
        w1 = uniform(low=0, high=1, size=1)
        u1 = uniform(low=0, high=1, size=1)
        Zi = jv(i, 2*r*np.sqrt(-v1 * np.log(w1))) * np.cos(i*theta+2*np.pi*u1)
        z_output.append(Zi)
    return np.sqrt(2) * np.sum(z_output)

Z = generate_z(example_1)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
    show=True, usecolormap=True,
    highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
                  width=1000, height=1000,
                  margin=dict(l=65, r=50, b=70, t=80), scene=dict(
        xaxis=dict(showticklabels=False),
        yaxis=dict(showticklabels=False),
        zaxis=dict(showticklabels=False),
    ))
fig.show()

(ii) Example 3.3 (continued), Case II: 0 < \nu < 1

We want \begin{equation} y_{1} = \Gamma(0.5, 1) \end{equation}

\begin{equation} y_{2} = U_{2^{0.5}, 8^{0.5}} \end{equation}\begin{equation} v_{1} = \frac{y_{1}^{0.5}}{y_{2}} \end{equation}\begin{equation} w_{1} = U_{0, 1} \end{equation}\begin{equation} u_{1} = U_{0, 1} \end{equation}

def example_2(i, j):
    r = R[i, j]
    theta = T[i, j]
    
    z_output = []
    for i in range(simulation_number):
        y1 = gamma(shape=0.5, scale=1, size=1)
        y2 = uniform(low=(2**0.5), high=(8**0.5), size=1)
        v1 = (y1 ** 0.5) / y2
        w1 = uniform(low=0, high=1, size=1)
        u1 = uniform(low=0, high=1, size=1)
        Zi = jv(i, 2*r*np.sqrt(-v1 * np.log(w1))) * np.cos(i*theta+2*np.pi*u1)
        z_output.append(Zi)
    return np.sqrt(2) * np.sum(z_output)

Z = generate_z(example_2)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
    show=True, usecolormap=True,
    highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
                  width=1000, height=1000,
                  margin=dict(l=50, r=40, b=80, t=75), scene=dict(
        xaxis=dict(showticklabels=False),
        yaxis=dict(showticklabels=False),
        zaxis=dict(showticklabels=False),
    ))
fig.show()

(iii) Example 3.3 (continued), Case III: \nu = 1 or \nu > 1

We want \begin{equation} y_{1} = \Gamma(3, 1) \end{equation}

\begin{equation} y_{2} = U_{\frac{1}{64}, \frac{1}{4}} \end{equation}\begin{equation} v_{1} = y_{1}^{0.5} * y_{2}^{0.25} \end{equation}\begin{equation} w_{1} = U_{0, 1} \end{equation}\begin{equation} u_{1} = U_{0, 1} \end{equation}

def example_3(i, j):
    r = R[i, j]
    theta = T[i, j]
    
    z_output = []
    for i in range(simulation_number):
        y1 = gamma(shape=3, scale=1, size=1)
        y2 = uniform(low=(1/64), high=(1/4), size=1)
        v1 = (y1 ** 0.5) * (y2 ** 0.25)
        w1 = uniform(low=0, high=1, size=1)
        u1 = uniform(low=0, high=1, size=1)
        Zi = jv(i, 2*r*np.sqrt(-v1 * np.log(w1))) * np.cos(i*theta+2*np.pi*u1)
        z_output.append(Zi)
    return np.sqrt(2) * np.sum(z_output)

Z = generate_z(example_3)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
    show=True, usecolormap=True,
    highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
                  width=1000, height=1000,
                  margin=dict(l=30, r=40, b=50, t=60), scene=dict(
        xaxis=dict(showticklabels=False),
        yaxis=dict(showticklabels=False),
        zaxis=dict(showticklabels=False),
    ))
fig.show()

(iv) Example 3.5 (continued)

def besselnum(n):
    results = []
    i = 0
    while i < n:
        u = uniform(0, 1)
        ug = uniform(0, 1)
        g = 2 * ug/(1-ug)
        if 3 * (0.5 / ((1+0.5*g) ** 2) * u) <= 2/g*(jv(g, 1) ** 2):
            i+=1
            results.append(g)
    return results

def example_4(i, j):
    r = R[i, j]
    theta = T[i, j]
    
    u = uniform(low=0, high=1, size=simulation_number)
    v = besselnum(simulation_number)

    z_list = []
    for i in range(simulation_number):
        Z = jv(i, r * v[i]) * np.cos(i*theta+2*np.pi*u[i])
        z_list.append(Z)

    return np.sqrt(2)*sum(z_list)

Z = generate_z(example_4)
fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
    show=True, usecolormap=True,
    highlightcolor="limegreen"))
fig.update_layout(title='', autosize=False,
                  width=1000, height=1000,
                  margin=dict(l=50, r=50, b=55, t=80), scene=dict(
        xaxis=dict(showticklabels=False),
        yaxis=dict(showticklabels=False),
        zaxis=dict(showticklabels=False),
    ))
fig.show()

(v) An example from Table 1

For simulation of

\begin{equation} Z(x) = \sqrt(2) \sum_{n=0}^{\infty}J_{n}(rV_{n})cos(n\pi + 2\pi U_{n}) \end{equation}

We being with simulating uniform random variables \begin{equation} U_{0}, U_{1}, U_{2} \sim [0,1]\end{equation}

rand_uniform = uniform(low=0, high=1, size=3) # simulate 3 random uniform
u0, u1, u2 = rand_uniform

The distribution function $ F(u) $ has to be chosen from say Table 1.

Consider the function on the second row of Table 1. When $d=2, \alpha=1$ we have \begin{equation} f(u) = \frac{1}{2} * u * exp(\frac{-u^2}{4}), u > 0 \end{equation} and \begin{equation} F(u) = 1 - exp(\frac{-u^2}{4}), u > 0 \end{equation}

Thus, we can use the inverse transformation method to generate random numbers from $F(u)$, where the inverse function of $F(u)$ is \begin{equation} F^{-1}(u) = 2 * \sqrt{-\log(1-u)} \end{equation}

For a uniform random variable $U$ on (0, 1), $U$ and $1-U$ have the same distribution. Thus, a code to generate one $F(u)$ random number is

\begin{equation}v = 2 * \sqrt{\log{U_{0}}} \end{equation}

rand_v_uniform = uniform(low=0, high=1, size=3)
v0, v1, v2 = 2 * np.sqrt(-np.log(rand_v_uniform))

Demonstrate Z(x) in formula 3.7 with an example, with n=2

np.sqrt(2) * (
    jv(0, r*v0) * np.cos(0*theta+2*np.pi*u0) + \
    jv(1, r*v1) * np.cos(1*theta+2*np.pi*u1) + \
    jv(2, r*v2) * np.cos(2*theta+2*np.pi*u2)
)

-0.05915042467644759

We show that the following formula is equivalent

def func_v(x):
    return 2 * np.sqrt(-1 * np.log(x))

np.sqrt(2) * sum(
    [
        jv(i, r*func_v(rand_v_uniform[i])) * \
        np.cos(i*theta+2*np.pi*rand_uniform[i])
        for i in range(0, 3, 1)
    ]
)

-0.05915042467644759

def generate_function(n, r, theta):
    n_random_uniform = uniform(size=(n+1))
    n_random_v_uniform = uniform(size=(n+1))
    
    return np.sqrt(2) * sum(
            [
                jv(i, r*func_v(n_random_v_uniform[i])) * \
                np.cos(i*theta+2*np.pi*n_random_uniform[i])
                for i in range(0, n+1, 1)
            ]
        )

# now calculate z for each r, theta

z_list = []
for theta in theta_list:
    for r in r_list:
        z_list.append(generate_function(200, r, theta))
        
## the output here is m_theta by m_r

Z = np.reshape(z_list, (m, m))

fig = go.Figure(data=[go.Surface(x=X, y=Y, z=Z)])
fig.update_traces(contours_z=dict(
    show=True, usecolormap=True,
    highlightcolor="limegreen",
    project_z=True))
fig.update_layout(title='Isotrophic Random Fields', autosize=False,
                  width=1000, height=1000,
                  margin=dict(l=60, r=50, b=60, t=90), scene=dict(
        xaxis=dict(showticklabels=False),
        yaxis=dict(showticklabels=False),
        zaxis=dict(showticklabels=False),
    ))
fig.show()